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Monte Carlo simulation techniques, like simulated annealing and parallel tempering, are often used to evaluate low- 
temperature properties and find ground states of disordered systems. Here we compare these methods using direct cal- 
culations of ground states for three-dimensional Ising diluted antiferromagnets in a field (DAFF) and three-dimensional 
Ising spin glasses (ISG). For the DAFF, we find that, with respect to obtaining ground states, parallel tempering 
is superior to simple Monte-Carlo and to simulated annealing. However, equilibration becomes more difficult with 
increasing magnitude of the externally applied field. For the ISG with bimodal couplings, which exhibits a high 
degeneracy, we conclude that finding true ground states is easy for small systems, as is already known. But finding 
each of the degenerate ground states with the same probability (or frequency), as required by Boltzmann statistics, 
is considerably harder and becomes almost impossible for larger systems. 
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1. Introduction 

One of the major current challenges in statistical physics is the study of random systems. Systems with 
quenched disorder like spin glasses and random field systems^ have attracted much attention during the 
past two decades. One is especially interested in numerically investigating the glassy low-temperature and 
ground-state behavior. In addition to the direct application, the calculation of ground states can be used to 
generate low-lying exited states^''^ or domain walls in systems to test whether the given system exhibits an 
ordered phase at non-zero temperature^'^'^. Recently, due to increased interactions between physicists and 
computer-scientists, many new ground-state techniques have been developed^ in this field. 

Which algorithm is suitable depends on the nature of the problem. In computer science there are two 
major types of optimization problems^: easy and hard ones. For easy problems there exist algorithms 
which can find an exact optimum while the running time of the algorithm increases only polynomially as a 
function of system size. For hard problems where no such algorithms are known so far, the running times 
of exact algorithms grow at least exponentially with the system size. In statistical physics, the ground 
state computation of three-dimensional random field Ising systems is easy, while obtaining ground state for 
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three-dimensional Ising spin glasses is hard. Both models are treated in this paper. 

For easy problems one can apply exact algorithms to study large systems. But for hard problems no 
fast algorithms arc available. Thus, very often heuristics and approximation methods arc applied. Many 
widely used techniques are based on Monte Carlo simulations^ '^'^'^^ (MC), which allow one to study a 
thermodynamical model at a given temperature T. To obtain low-temperature and ground states, many 
variants have been developed, such as simulated annealing^^ (SimA) which mimics slowly cooling the system 
and parallel tempering^'^'^'^ (PT), where a system is simulated at several temperatures in parallel. 

In this paper we focus on parallel tempering, a method that has been recently applied widely to spin 
glasses^^'^^'^^, random field systems^^, clusters/molecules^^ and various other systems^°. To our knowledge, 
so far only comparisons of different MC techniques with PT have been performed. This has the disadvantage 
that for all of these techniques, including PT, the proof of equilibration depends on the model and is difficult 
to obtain. Here we assess how powerful PT is in generating ground states. We apply PT to the diluted 
(Ising) antiferromagnet in a field (DAFF) and to the Ising spin glass (ISG). Two questions are asked: 

• How often are ground states obtained ? 

• Is the ground-state sampling statistically correct in degenerate systems? 

For the DAFF, there exists a fast exact ground-state algorithm, hence we can study whether PT is able to 
obtain exact ground states and compare to other methods. This kind of comparison has the advantage that 

the correct result is known exactly and there are no problems with equilibration for the exact ground-state 
calculations. This is different to previous work, where only different approximation techniques have been 
compared. We demonstrate here that PT is able to generate exact ground states for the DAFF, and that 
is it superior to simple MC and to SimA. Still, the method has its limitations: with increasing size of the 
external field, obtaining ground states for the DAFF becomes increasingly difficult. 

We report another restriction of the method in the second part, where PT is applied to the three- 
dimensional ISG with bimodal distribution of the interactions. This system has a large ground-state degen- 
eracy. For an unbiased sampling of the ground states, each ground state must be present in any result with 
the same probability. By studying small systems that have only few different ground states and where one 
can do many independent runs, we show that PT, while able to obtain true ground states for the ISG, has 
a harder time to get an unbiased sampling. Hence, for larger systems where the degeneracy is larger and 
the algorithm is subsequently slower, additional measures must be undertaken^^ to have an unbiased result. 
This finding also applies to studying the low-temperature (i.e. non-zero) behavior in general. 

The rest of the paper is organized as follows. First we describe the models that are studied and the 
observables investigated in this paper. Next, the different algorithms used arc explained. In the two following 
sections, the results and comparisons for the application of the algorithms to the DAFF and the ISG are 
presented. In the last section we conclude with a summary. 

2. Models and Observables 

To assess the efficiency of PT, SimA and simple Monte Carlo, we apply them to the DAFF and the ISG. In 
this section both models and background information are presented. 

An example for a physical realization of the DAFF is the crystal FexZni_xF2. Besides growing prac- 
tical applications^^ of this system, many fundamental questions arc still unsolved. For instance, at the 
antiferromagnet-paramagnet transition some experiments^^'^^ find a logarithmic divergence of the specific 
heat, corresponding to a specific heat exponent a = 0. However, numerical simulations^^'^^'^'' give a negative 
value of a « —0.5. Also a non-diverging susceptibility was found^* using exact ground-state calculations. 
Therefore the nature of the transition remains unanswered. 

Two theoretical models have been proposed: the random field Ising model (RFIM) and a direct model 
for the DAFF. Both models are expected to be in the same universality class^^'^°, hence results from both 
should be equivalent. Since the RFIM is easier to simulate than the DAFF, it has been studied more despite 
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the fact that the DAFF is closer to the physical system. In addition, disagreements about the universality 
have been found^^'^^'^^ at T = where the critical exponents of the ±A RFIM and the DAFF seem to be 
significantly different. Furthermore, there is still no finite-temperature study for large system sizes. This 
was the initial motivation to start the current work. 
The DAFF is given by the following Hamiltonian: 

H = -j'^eiCjSiSj - H'^CiSi , (1) 

ihj) ' 

where the sum is over the nearest neighbors on a cubic lattice of size Lx Lx L. The 5j = ±1 represent Ising 
spins. Periodic boundary conditions arc applied. As we want antifcrromagnctically coupled spins, we choose 
J = — 1 and H represent the externally applied magnetic field. The dilution is introduced by randomly 
setting ei = for a vacant site and otherwise = 1. A realization of the disorder is given by a set {ci} 
of fixed (i.e. quenched) values. In this work we use a dilution of 50%, i.e., we have the same number of 
randomly chosen empty and filled sites, to ensure that both occupied and vacant sites are percolating. For 
the final results, an average over the disorder has to be performed. 

For the DAFF, we are studying the physical quantities given below. They are obtained from the recorded 
values of the energy per spin and the staggered magnetization per spin M . The energy per spin is given by 
the thermal and disorder average of the Hamiltonian of the system E{S)/N, where A'^ is the number of spins 
for 50% dilution). The microscopic staggered magnetization is given by 

1 ^ 

M = -Y.^-ir^'+^Si. (2) 

Here, X, y, and z are the spatial coordinates of the spin 5*^. The factor (— l)^+f+^ takes the order of two 
sub-lattices into account and ensures that the above magnetization is unity when the system is ordered. 
We evaluate the average staggered magnetization defined as 

Hav = [\{M)\U ■ (3) 

Here (. . .) represent a thermal average whereas [. . .]av represents an average over disorder. 

The staggered magnetization is the order parameter for the DAFF. The canonical phase diagram of the 
DAFF in higher than two dimensions is shown in Fig. 1. For low temperatures and small randomness, the 
antifcrromagnctic couplings dominate, hence the system exhibits a long-range order. For higher temperatures 
T > Tc{H), where the entropy dominates, or for higher fields, where the spins are predominately aligned 
parallel to the field, the system is paramagnetic. 

In addition, we evaluate the average specific heat: 

[CU=N{[{E')U-[{E)X^)/T^ (4) 

the susceptibility 

[xU = N{[{M')U-[{\M\fU), (5) 

and the Binder cumulant^^, which related to the ratio between the fourth moment of the magnetization 
divided by the second moment squared, i.e. 



av 



Furthermore, we study the Edwards- Anderson (EA) Ising spin glass. The Hamiltonian of the EA model 

is given by 

W = - ^ JrjSiSj , (7) 

(hi) 
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Figure 1: A sketch of the phase boundary of the DAFF. The antiferromagnetic phase is denoted by "AF" 
and the paramagnetic phase by "P" . 

where the Ising spins 5, = ±1 lie on a hypercubic lattice in three dimensions with N = sites, and Jij are 
nearest-neighbor interactions chosen according a given probability distribution. We study the Gaussian case 
as it is known to have a unique ground state (up to a global symmetry) as well as the bimodal case which 

is highly degenerate* 

For the ISG, we study the ground state energy and record histograms of how often each of the degenerate 
ground states is found by performing several simulations with the same disorder realization but different 
initial spin configurations. 

3. Algorithms 

We study the DAFF and the ISG with 5 different algorithms. Three of them are simulation methods at 
finite temperature: simple Monte Carlo, simulated annealing and parallel tempering Monte Carlo. For the 

ground state calculations, wc use exact graph-theory based techniques for the DAFF and a heuristic method, 
the genetic cluster-exact approximation algorithm, for the ISG. 

A review on recent simulation techniques at finite temperatures can be found in Ref. Here we 
just state the fundamental notions behind the applied techniques. The basic idea of the Monte Carlo 
approach is that a random walk in configuration space is performed, leading to a stationary distribution of 
the configurations S which is the Boltzmann distribution at fixed temperature T. At lower temperatures, the 
algorithm freezes, hence calculating ground states is difficult, especially for random systems like spin glasses 
and diluted antiferromagnets. 

A more sophisticated approach to obtain low-temperature states is the simulated annealing (SimA) 
technique^^. The basic idea is to mimic the experimental techniques that are used to determine the low- 
temperature behavior of a probe: the system is equilibrated at high temperature, well above the phase 
where the glassy behavior occurs. Then the system is slowly cooled down to the desired temperature. But, 
it remains difficult to obtain equilibration. The reason is that one might still be caught in local minima. 
Since the temperature is monotonically decreased, it is hard to escape a minimum once the temperature is 
already low. 

For this reason another method has been invented, the parallel tempering (PT) approach^'^''^'*. The 
basic idea is that the temperature is not only reduced, but it is also allowed to raise from time to time to 
escape local minima of the free energy landscape. This algorithm is realized by keeping Nt different and 
independent configurations at temperatures Ti„i„ = Ti < T2 < . . . , < Tjv^ = Tmax For each configuration, 

*Spin glasses have been studied widely in statistical physics, for more information we refer to the literature^*'^. 
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several MCS are performed. Then the configurations are aUowed to exchange their positions in the fist, i.e. 
to exchange their temperatures. The swapping probability is chosen, such that for joint distribution over all 
temperatures the joint Boltzmann distribution is obtained. 

To assess the ability of PT and SimA to generate ground states, we compare with exact ground states. 
The ground state calculation of the DAFF can be solved in a time that is polynomial in the number of spins. 
To calculate the exact ground states at given randomness {e,} and field H, algorithms^^'^^'^®''' from graph 
thcory^°'^^'^^ are applied. The calculation works by transforming the system into a network''^, calculating 
the maximum flow in polynomial time^^'^^'^^'^'''^^, and finally obtaining the spin configuration {Si] from the 
values of the maximum flow in the network. The running time of the latest max;imum-flow method has a 
peak near the phase transition and diverges^^'^" there like 0{L'^^^), hence very large systems like N = 100^ 
can be studied easily. Here, where we are interested in comparison with results from Monte Carlo techniques, 
we restrict ourself to small systems up to = SG'^. 

Finally, for two-dimensional spin glasses without external field with periodic boundary conditions at 
most in one direction, efficient polynomial algorithms for the calciilation of exact ground states arc available. 
Recently, results for systems of size 1800 x 1800 have been obtained^^. On the other hand, the calculation of 
ground states for three-dimensional ISGs belongs to the class of the NP-hard problems^^, i.e., only algorithms 
with exponentially increasing running time are available. Thus, only small systems can be treated^^. The 
basic method used here is a combination of a genetic algorithm and the cluster-exact approximation (CEA) 
technique^'* which allows to obatin ground states up to size N = 14'^. A pedagogical presentation of this 
combination of methods can be found elsewhere ^. 

4. Results for the DAFF 

We have studied the DAFF with system size L = 12. Our initial motivation has been to use PT to go to 
much larger systems sizes. Simulations for sizes L — 24, and L = 36 have been performed, but only for 
H = 0.4. Our experience is that equilibration is too diflicult, i.e., it is impossible to obtain reliable results 
for the larger sizes. Since we arc only interested in comparing different algorithms, small sizes arc sufficient. 

For parallel tempering to work well one has to choose carefully a set of temperatures to ensure that the 
acceptance ratios for the global moves are constant (as a function of T), preferentially near 0.5. The number 
of temperatures depends on the size of the lattice and the lowest temperature studied. We have used the 
method introduced by Hukushima^^ where one recursively adapts the temperature set until the acceptance 
ratios are about 0.5 for all temperature pairs. 

The resulting number Nt of temperatures, the minimum and maximum temperatures, the number of 
samples (i.e. disorder realizations) A'^gamp, the total number of sweeps Agwccp are shown in Table . 

The simulations are started with a random spin configuration. We have checked that the acceptance 
ratios are almost fiat [see Fig. 2(a)]. In Fig. 2(b) we show data for different observables for the DAFF as 
a function of Monte Carlo steps in order to see that the system properly equilibrates when using PT. The 
equilibration time N^q is the number of MCS needed at the lowest temperature such that all observables 
become time-independent, when averaged over the last half of a MC run. The numbers for Aeq given in 
Table are only an upper bound with an accuracy of a factor two. 

All these parameters are the same for the three methods used in order to compare the efficiency of the 
algorithms. This means for simple MC wc have simiilatcd the system for all temperatures Ti independently. 
For SimA, we have started to simulate the system with temperature Iat^ and cooled stepwise down to Ti. 
Per temperature Ngweeps MCS were performed. This means we have not used the usual exponential cooling 
schedule. Please note that the Hukushima scheme results in more temperatures in the low temperature 
regime, hence the difference to the exponential scheme is not too large. 

In Fig. 3 the results at H = 0.4 for the specific heat, the susceptibility, the staggered magnetization and 
for the Binder parameter are shown as a function of the temperature T for all three methods used. The 
average is taken for MC times larger than A^eq and over all realizations of the disorder. The error bars are 
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Parameters of the Simulation 



H 


Nt 




^ samp 


^ sweep 




0.4 


40 


[0.0128, 4.2613] 


260 


524288 


65536 


0.8 


42 


[0.9570, 4.2613] 


150 


524288 


65536 


1.2 


47 


[0.0328, 3.7144] 


120 


524288 


65536 


1.6 


50 


[0.0497, 4.0419] 


140 


131072 


32768 


2.0 


46 


[0.0003, 4.0419] 


400 


65536 


16384 



Table 1. Parameters of the simulation for the DAFF: H is the externally applied magnetic field, and Nt the number of 
temperatures used in the parallel tempering method. A^aamp is the total number of disorder realizations, Nnwcap denotes 
the number of MCS per sample and temperature, while A'cq is the equilibration time after which the measurements 
have started. 




Fig. 2. Panel (a) shows the acceptance ratios that are obtained after the simulation for the DAFF. Note that they 
are around ~ 0.5 and almost temperature independent. 

In panel (b) we show different observables as a function of MCS for L = 12 (DAFF). One can clearly see that the 
staggered magnetization [m]av, specific heat [Cjav, absolute value of the energy per spin |[e]ac| and susceptibility [x]av 
become independent of the number of MCS at the same equilibration time. The data is for H = 0.4 and T = 1.81. 
This temperature, close to the critical, shows the worst equilibration. The data for the susceptibility have been divided 
by 20 for better viewing. 

from sample to sample fluctuations. One can clearly observe that for large temperatures all three methods 
give the same result, while for low temperatures the simple IVIctropolis approach exhibits large unphysical 
deviations, especially for the magnetization. Hence, for the rest of the chapter we focus on the comparison 
between PT and SimA. 

In Fig. 4 the staggered magnetization obtained from these two methods is compared for larger values of 
the field H = 0.8 and H = 1.2, respectively. With growing strength of the field, the deviations between both 
methods increase. 

We want to understand the cause of the difference. For this reason, we apply the exact graph-theoretical 
algorithm to obtain exact ground states and compare to the results obtained by PT and SimA. In Fig. 5 
the averages /i and /2 are presented, /i is the sample average of an indicator function, which is 1 if for a 
sample at least once a ground state was found at a given temperature, and if a ground state has never been 
visited. /2 is the average fraction of the simulation time a system has been in a ground state at temperature 
T. Note that both quantities must converge to 1 when T ^ 0, and that fi{T) > f2{T) always holds. 

It is clearly visible that SimA fails in finding the true ground state for a large fraction of all realizations. 
Furthermore, both /i and /2 seem not to converge to 1 for T ^ 0. With increasing size of the external field 
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Fig. 4. Zoom on the results for the staggered magnetization at low temperatures for L = 12 with the applied fields 
H = 0.8 (a) and H = 1.2 (b) for the DAFF. Although for the rest of the range the result is almost the same in both 
methods for low temperatures, the results from parallel tempering are more reliable. 
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Fig. 5. This figure is a comparison between parallel tempering (PT) and simulated annealing (SimA) for various 
equilibration criteria for the DAFF; the average over of the disorder of the times the simulation finds the ground state 
at least once in the range of measurements (/i), and the average over the disorder and over the range of measurements 
of the times the system is in a ground state (/2)- 
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the ability of SimA to find ground states decreases, except for H = 2.0 where the system exhibits a large 
ground-state degeneracy, that seems to facilitate the ground-state search. Also, the results of SimA for /i 
exhibits strong fluctuations, another indication how unreliable the result is. 

On the contrary, the data for PT behaves as expected, both /i and /2 converge to 1 for T ^ 0. This 
in turn shows that it is very likely that indeed PT has obtained the true thermodynamic average, hence 
it is superior to SimA (and simple MC). Again, finding ground states becomes more difficult with growing 
size of the field, except the highly degenerate case H = 2.0. A possible explanation for this fact can be the 
existence of temperature "chaos" . This means that for neighboring temperatures, the free-energy landscape 
looks very different. Because PT relies on the fact that neighboring temperatures yield a similar behavior, 
the existence of chaos would prevent a quick equilibration. So far mostly spin glasses have been investigated 
regarding the existence of chaos. For zero field no definite answer on the existence of chaos has been given, 
but chaos seems to be weak and/or difficult to observe 55,56,57,58^ p^^. ^j-^^ RFIM, which is very similar to 
the DAFF, chaos has been studied only for small random perturbations of the bonds. In this case, for the 
three-dimensional model, chaos seems to be marginal^^. For the case of temperature chaos in the DAFF, we 
believe that due to the additional field term in the Hamiltonian, the effect might be much larger for systems 
in an external field. Our result that finding ground states becomes increasingly difficult with H > indicates 
that this might be indeed the case. 

Note that using this comparison of the ability to obtain ground states, the deviations between the different 
algorithms are clearer than by simply comparing physical quantities. Especially for low values of the applied 
field H, PT and SimA give comparable results (cf. Fig. 3). The reason is that for low values of H, the 
physical quantities for ground states and excited states do not difi'er to much, hence the failure of SimA is 
not revealed. 

In this section we have demonstrated that parallel tempering is superior to simple MC simulations and 
simulated annealing, for moderate system sizes. The reason being that PT is much better in sampling the 
low-temperature region, especially the ground states, as demonstrated by comparison with exact ground 
states obtained by a graph theoretical algorithm. In the next section we compare PT and SimA for another 
model where no fast exact ground-state algorithms are available. 

5. Ground-state Statistics for the 3D Edwards- Anderson Ising Spin Glass 

In this brief section, we apply PT and SimA to the three-dimensional Edwards- Anderson (EA) Ising spin glass 
with a Gaussian respectively bimodal distributions of the interactions. We show that PT is able to calculate 
true ground states and that we obtain the statistics of the ground states for the degenerate bimodal systems. 
For the results given below, we have also checked with the genetic cluster-exact approximation method 
explained above that indeed ground states are found. 

To determine the energy of the system for a given set of bonds, we equilibrate for a given number of 
MCS, in this case 10^ per temperature set, and take a snapshot of the system configuration. Figures 6 (a) 
- (c) show data for the probability to find a given energy P{E) at temperature T » O.lTc (Tc « 0.95)^° of 
three single realizations with Gaussian distributed bonds with L = 4, 6, and 8. We choose the Gaussian 
case first as we have a strong test^^ to see when the system has equilibrated. We see that for low enough 
temperatures, in this particular case, the ground state Eq is the most probable state for the system. For the 
simulations we have used a set of 18 temperatures in a range [O.I, 2]. 

We have again measured the equilibration time^^ iVeq i.e. the number of MCS needed such that all 
measured quantities become independent of the simulation time. Interestingly, it increases for this small 
system sizes with a power law as a function of the number of spins, as can be seen in Fig. 6 (d). Since the 
spin-glass ground-state calculation is computationally a hard problem^^, we expect that for larger system 
sizes, the running time increases exponentially with system size. This has been observed for algorithms that 
are especially designed to compute spin glass ground states^^'^^. Note that PT is slower than those special 
algorithms, but on the other hand is applicable to a variety of models due to its general structure and very 
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simple implementation. 



~l ' I ' I ' 
L = 4 Gaussian (a) 
Eq -1.712031 
T = 0.1 
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-1.767 
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Eq = -1.705370 
T = 0.1 
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Fig. 6. Panels (a) — (c) show data for the probability to find a given energy P{E) for the EA model with Gaussian 
bonds for three single realizations of size L = 4, 6, and 8. One can see that the ground state energy Eq is the most 
probable state for T = 0.1. Panel (d) shows the total number of MCS needed to (including the number of temperatures 
used and spins in the system) evaluate one ground-state realization as a function of the number of spins N. A clear 
power-law increase with a power of 2.61 can be seen. 



To further illustrate the use of PT and SimA to determine ground states for spin-glass systems, we have 
applied them to the highly degenerate bimodal ISG. For the Gaussian case above, we have equilibrated the 
system at low temperatures, then we have checked that indeed ground states are found. Now we proceed in a 
different way. We make the simulations for different number of MCS and investigate directly whether ground 
states are obtained by comparing with the results from the genetic CEA approach. If indeed ground states 
are found, we furthermore check if finding ground states automatically means that the system is equilibrated. 
This is tested by measuring the frequency different degenerate ground states are found. 

In this case we have used a set of Nt = 10 temperatures in the range [0.35, 2]. We present in Fig. 7 (a) 
the ground-state statistics for a specific choice of interaction bonds for the bimodal case calculated with PT 



for TV,, 



10^ MCS per replica, and iVn 



10 independent runs. The system we have studied has 64 



distinct ground states. For a thermodynamically correct sampling, each ground state has to be found with 
the same frequency^^. Hence, equilibration can be tested unambiguously by counting how often each ground 
state is found. As we can see, the distribution of ground states is not completely flat, although the algorithm 
often finds a ground state, hence is already very efficient. At least the results are better than those obtained 
by means of simulated annealing [see Fig. 7 (b)] for the same total number of MCS (10'^) and within the 
same range of temperatures, but with an exponential temperature schedule T„+i ~ aT„ 
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Fig. 7. Panel (a) shows the number of hits for a specific ground-state configuration for the bimodal EA spin glass 
for T = 0.35 using parallel tempering (PT) Monte Carlo. Data for A^sweep = 10^ MCS, Nt = 10 temperatures and 
Nruns = 10^ independent runs. For our particular choice of bonds, we find 64 distinct ground-state configurations. 
Note that the distribution of states, while not completely flat, gives almost equiprobable ground-state conflgurations 
with Q = 0.09. This is in contrast to panel (b) where we show data obtained from simulated annealing (SimA) where 
Q = 0.18 for a cooling rate a = 0.9827, which corresponds to 10^ MCS. 




20 40 60 20 40 60 



GS GS 

Fig. 8. Same as Fig. 7 but for 10* MCS. One can see that for parallel tempering (PT) Q = 0.04 whereas for simulated 
annealing (SimA) Q = 0.08 for a cooling rate a = 0.9983 which corresponds to 10* MCS. 
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In Fig. 8 we show the same results as in Fig. 7 but for 10^ MCS. One can clearly see that fluctuations 
are reduced. 



Ground-state statistics for T = 0.35 



-^sweep 


Nt 


N{Eo) 


e{GS) 




Q 


10^ 


10 


68987 


1061.34 


490.04 


0.46 


10^ 


10 


99973 


1538.05 


133.91 


0.09 


lO'^ 


10 


99972 


1538.03 


66.22 


0.04 


10* 


10 


99970 


1538.00 


63.78 


0.04 


lO'' 


10 


99976 


1538.09 


66.87 


0.04 



Table 2. Ground-state statistics for T = 0.35 in the bimodal EA model. For increasing number of MCS {Nsy,eep), 
the number of times we reach the ground state N{Eo) increases. Also shown are the mean number of times a specific 
ground state configuration is reached [e(GS) = N{E(j)/6'l], as well as its standard deviation [(t(GS')] and relative 
fluctuations [Q = a(GS)/e(GS)]. Note the number of MCS quoted is jjer replica. Hence, the total number of MCS of 
the simulation is given by N'y ■ Ns-weep, where is the number of temperatures (replicas) used. Data for A^samp = 10^ 
disorder realizations. 



Ground-state statistics for T = 0.65 



^ sweep 






e{GS) 


aiGS) 


Q 


10^ 


10 


10322 


158.80 


70.27 


0.44 


10^ 


10 


87951 


1353.09 


103.92 


0.08 


lO'' 


10 


88446 


1360.71 


50.47 


0.04 




10 


88298 


1358.43 


59.10 


0.04 


10*= 


10 


88673 


1364.20 


57.48 


0.04 



Table 3. Same as Table 2 but for T = 0.65. 



In Tables 2 and 3 we show the dependence of the number of times we reach the ground state as a function 
of MCS for T = 0.35 and 0.65, respectively. (For the bimodal case w 1.15)^^''^''. While the system is 
known to be properly equilibrated for A^swccp = 10'"' MCS (per temperature set) for T = 0.35, we see that 
for shorter equilibration times we still reach the ground state several times. We also calculate the mean 
number of times a ground state is reached [e{GS)] and estimate a relative fluctuation Q by dividing the 
standard deviation [(t(G'S')] by the mean. The distribution of ground states is fully equiprobable in the limit 
Q — *■ 0.025. This can be understood as follows: drawing from G numbers randomly (in this particular case 
we have G = 64 ground states) with all being equiprobable with probability P = 1/G follows a binomial 
distribution. If we repeat this times we obtain a standard deviation a = ■\/Np{l — p) and expectation 
value e = Np. This means that Q = \/ {p — l)/{Np), and for our case with G = 64 and N = 10^ (= A?runs) 
we obtain Qth = 0.025. 

As we can see, fluctuations decrease dramatically for both temperatures when the number of MCS is 
increased but saturates at Q = 0.04. We see that we are rather close; to the theoretical value showing 
that while the fluctuations are very small it would require prohibitively longer equilibration times to flnd 
Qth = 0.025. 

Here again, parallel tempering shows an improvement over simulated annealing when trying to determine 

ground states and ground-state statistics. In particular, we sec that for small systems with moderate effort 
we can determine the ground state of a spin system with equiprobable ground-state statistics. The eflfort 
increases considerable for larger systems. There, not only the computational eff'ort to find some ground 
states is already significantly larger, but also the ground-state degeneracy grows exponentially with system 
size. For example, a system of size N = 6'^ can have 10^^ ground states, and it is impossible to generate a 
histogram as done for small system sizes. Hence, when evaluating ground-state properties, one cannot take 
the results from MC-like simulations directly. Finding ground states in equilibrium is more diSicult than 
just finding ground states. Therefore, one has to use long simulation times at low temperatures and wait for 
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measurable quantities to become time-independent, or additional algorithms'^ must be applied to guarantee 
a correct sampling. 

6. Summary and discussion 

In conclusion, we have studied the three-dimensional diluted Ising antiferromagnet in a magnetic field 
and the Edwards-Anderson Ising spin-glass using parallel tempering Monte Carlo, simulated annealing and 
simple Metropolis Monte Carlo. We have used the generation of ground states as an assessment tool to 
compare these three techniques. We find the following: 

• PT is superior to simulated annealing and simple Monte Carlo: More ground states are generated at 
very low temperature, where the system is preferentially in a true ground state for small systems. 

• For systems in an external field, generating ground states becomes more difficult, except if the degen- 
eracy is very large. We have also performed similar tests for the ISG where we have observed the same 
behavior (not shown here). A possible explanation could be the existence of stronger temperature 
"chaos" in comparison to the zero field case. 

• Although PT is superior to the other methods studied, generating ground states with the correct 
statistics is more difficult than simply generating ground states at all. For larger than tiny systems, 
it is impossible to assess directly whether the sampling is correct. Either additional methods must be 
applied or very long MC simulation times must be accepted. 

We conclude: for system where exact algorithms are available that run in polynomial time, like for 
the DAFF, it is recommended not to use a simple Monte Carlo approach (like PT) for the ground state 
calculation. The reason is that one can treat much larger system sizes by using the exact techniques than 
it is possible with MC, while MC does not guarantee true ground states. If one is interested not only in the 
exact zero-temperature behavior but in the non-zero low-temperature behavior, one should use a combination 
of the exact techniques and PT. When only applying Monte-Carlo techniques like PT, it seems unlikely that 
for the DAFF it is possible to equilibrate large systems (such as 100^ spins). In this aspect the DAFF is of 
similar complexity as spin glasses. This seems strange, because for the ground-state problem, the DAFF is 
much easier to simulated than ISGs. 

For the ISG and other problems, where no fast exact algorithms are available, currently PT is the best 
choice for an algorithm that is general purpose and simple to implement. But, as we have seen, one cannot 
say that PT is orders of magnitudes faster than standard algorithms used. It suffers from the same problems 
as SimA or simple MC, but equilibrates faster. Care has to be exercised before results are believed. As we 
have seen it is still easily possible that the true equilibrium behavior is not obtained. Several independent 
equilibration tests should always be applied. 
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